module ModuleDistributions_HSV

    use Globals
    use Functions_wa_HSV
    use Toolbox
    
contains
    
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! Distribution of tax-and-transfer rates, income, net worth
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    subroutine income_distribution

        implicit none

        ! Local variable declarations
        real(8), dimension(NS) :: yk, yl, ytot, tax_paid                                        ! capital, labor, and total income
        real(8), dimension(NS) :: ytot_sorted, yl_sorted, tax_paid_sorted, transfers_sorted     ! total income, labor income, tax paid, transfers (all sorted)
        real(8), dimension(NS) :: mu_sorted, mu_cumu, a_sorted, total_tax, total_tax_sorted     ! measure sorted, cum. measure, assets sorted, cap.+lab. tax, also sorted
        real(8), dimension(NS) :: av_tax_sorted, marg_tax_sorted, av_tax                        ! avg. tax sorted, marg. tax sorted, avg. tax 
        real(8), dimension(NS) :: marg_tax, dmarg_tax, dmarg_tax_sorted                         ! marg. tax, change in marg. tax, change in marg. tax sorted
        integer :: ind_sort(NS)                                                                 ! index for sorting
        integer :: Nstate                                                                       ! auxiliary variable equal to NS
        real(8), dimension(NS) :: muu, S1u, S2u                                                 ! auxiliary variables equal to measure, asset states, productivity states
        real(8) :: avg_tax_q(5), avg_transfers_q(5), inc_share(7), a_share(7)                   ! avg. tax-transfer by quintile, avg. transfer by quintile, inc. share by quintile and top 10 and top 1, asset shares
        real(8) :: avg_income_tax_q(5), a_share_top10, a_share_top1                             ! avg. lab. tax, asset shares top 10 and top 1
        real(8) :: avg_cap_income_tax_q(5), avg_total_tax_q(5), avg_tax_ind_q(5)                ! avg. cap. tax, avg. total tax, avg. tax other concept
        real(8) :: marg_tax_ind_q(5), dmarg_tax_ind_q(5)                                        ! marg. tax by quintile, change in marg. tax by quintile
        real(8) :: inc_share_top10, inc_share_top1                                              ! inc. share top 10 and top 1
        integer :: p_20_ind, p_40_ind, p_60_ind, p_80_ind, p_90_ind, p_99_ind, p_50_ind         ! indices for percentiles
        integer :: p_1_ind, p_5_ind, p_10_ind                                                   ! indices for percentiles
        integer :: iter, is                                                                     ! counter
        real(8) :: chit                                                                         ! auxiliary variable transfer
        real(8), dimension(NS) :: et, transfers                                                 ! auxiliary variable transfer, transer
        real(8), dimension(NS) :: income_tax, income_tax_sorted                                 ! income tax, and sorted
        real(8), dimension(NS) :: capital_income_tax, capital_income_tax_sorted, yk_sorted      ! cap. tax, sorted, cap. income sorted

        ! auxiliary variables
        Nstate = NS     ! number of individual states
        muu = mu        ! measure over individual states
        S1u = S(:,1)    ! assets for individual states
        S2u = S(:,2)    ! productivities for individual states

        ! Income variables
        yk = r*S1u          ! capital income
        yl = wge*S2u*h_pol  ! labor income
        ytot = yk+yl        ! total income
    
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        ! Tax and transfer rates
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

        ! Total taxes and transfers
        tax_paid = TAX(yl,yk)   ! total tax-and-transfer paid
        av_tax = tax_paid/ytot  ! total tax-and-transfer average rate

        do is = 1, Nstate
            marg_tax(is)  = dTAX_sc(yl(is),yk(is))      ! total tax-and-transfer marginal rate
            dmarg_tax(is) = ddTAX_sc(yl(is),yk(is))     ! total tax-and-transfer change in marginal rate
        enddo

        ! sort total income array
        ytot_sorted = ytot
        ind_sort = 0
        call hpsort_eps_epw(Nstate, ytot_sorted, ind_sort, 1D-16)
        ! sort measure and taxes paid arrays correspondingly
        mu_sorted = muu(ind_sort)                                   ! measure sorted by total income
        tax_paid_sorted = tax_paid(ind_sort)                        ! total tax-and-transfer sorted by total income
        av_tax_sorted = av_tax(ind_sort)                            ! average tax-and-transfer rate sorted by total income
        marg_tax_sorted = marg_tax(ind_sort)                        ! marginal tax-and-transfer rate sorted by total income
        dmarg_tax_sorted = dmarg_tax(ind_sort)                      ! change in marginal tax-and-transfer rate sorted by total income

        ! compute cumulative measure
        mu_cumu = 0d0
        mu_cumu(1) = mu_sorted(1)
        do iter = 2, Nstate
            mu_cumu(iter) = mu_cumu(iter-1) + mu_sorted(iter)
        end do

        ! Quintiles
        p_20_ind = minloc(abs(mu_cumu-0.20d0), dim=1)
        p_40_ind = minloc(abs(mu_cumu-0.40d0), dim=1)
        p_60_ind = minloc(abs(mu_cumu-0.60d0), dim=1)
        p_80_ind = minloc(abs(mu_cumu-0.80d0), dim=1)
        !print *, 'Top quintile cutoff = ', ytot_sorted(p_80_ind)/ymean

        ! compute average tax rates by quintile (CBO concept: total tax of the quintile divided by total income of the quintile)
        avg_tax_q(1) = sum(tax_paid_sorted(1:p_20_ind)*mu_sorted(1:p_20_ind))/sum(ytot_sorted(1:p_20_ind)*mu_sorted(1:p_20_ind))
        avg_tax_q(2) = sum(tax_paid_sorted(p_20_ind+1:p_40_ind)*mu_sorted(p_20_ind+1:p_40_ind))/sum(ytot_sorted(p_20_ind+1:p_40_ind)*mu_sorted(p_20_ind+1:p_40_ind))
        avg_tax_q(3) = sum(tax_paid_sorted(p_40_ind+1:p_60_ind)*mu_sorted(p_40_ind+1:p_60_ind))/sum(ytot_sorted(p_40_ind+1:p_60_ind)*mu_sorted(p_40_ind+1:p_60_ind))
        avg_tax_q(4) = sum(tax_paid_sorted(p_60_ind+1:p_80_ind)*mu_sorted(p_60_ind+1:p_80_ind))/sum(ytot_sorted(p_60_ind+1:p_80_ind)*mu_sorted(p_60_ind+1:p_80_ind))
        avg_tax_q(5) = sum(tax_paid_sorted(p_80_ind+1:Nstate)*mu_sorted(p_80_ind+1:Nstate))/sum(ytot_sorted(p_80_ind+1:Nstate)*mu_sorted(p_80_ind+1:Nstate))

        ! compute average tax rates, at the indiv level (mean of average tax-and-transfer rate)
        avg_tax_ind_q(1) = sum(av_tax_sorted(1:p_20_ind)*mu_sorted(1:p_20_ind))/sum(mu_sorted(1:p_20_ind))
        avg_tax_ind_q(2) = sum(av_tax_sorted(p_20_ind+1:p_40_ind)*mu_sorted(p_20_ind+1:p_40_ind))/sum(mu_sorted(p_20_ind+1:p_40_ind))
        avg_tax_ind_q(3) = sum(av_tax_sorted(p_40_ind+1:p_60_ind)*mu_sorted(p_40_ind+1:p_60_ind))/sum(mu_sorted(p_40_ind+1:p_60_ind))
        avg_tax_ind_q(4) = sum(av_tax_sorted(p_60_ind+1:p_80_ind)*mu_sorted(p_60_ind+1:p_80_ind))/sum(mu_sorted(p_60_ind+1:p_80_ind))
        avg_tax_ind_q(5) = sum(av_tax_sorted(p_80_ind+1:Nstate)*mu_sorted(p_80_ind+1:Nstate))/sum(mu_sorted(p_80_ind+1:Nstate))

        ! compute marginal tax rates by quintile (mean of marginal tax-and-transfer rate)
        marg_tax_ind_q(1) = sum(marg_tax_sorted(1:p_20_ind)*mu_sorted(1:p_20_ind))/sum(mu_sorted(1:p_20_ind))
        marg_tax_ind_q(2) = sum(marg_tax_sorted(p_20_ind+1:p_40_ind)*mu_sorted(p_20_ind+1:p_40_ind))/sum(mu_sorted(p_20_ind+1:p_40_ind))
        marg_tax_ind_q(3) = sum(marg_tax_sorted(p_40_ind+1:p_60_ind)*mu_sorted(p_40_ind+1:p_60_ind))/sum(mu_sorted(p_40_ind+1:p_60_ind))
        marg_tax_ind_q(4) = sum(marg_tax_sorted(p_60_ind+1:p_80_ind)*mu_sorted(p_60_ind+1:p_80_ind))/sum(mu_sorted(p_60_ind+1:p_80_ind))
        marg_tax_ind_q(5) = sum(marg_tax_sorted(p_80_ind+1:Nstate)*mu_sorted(p_80_ind+1:Nstate))/sum(mu_sorted(p_80_ind+1:Nstate))

        ! compute change in marginal tax rates by quintile
        dmarg_tax_ind_q(1) = sum(dmarg_tax_sorted(1:p_20_ind)*mu_sorted(1:p_20_ind))/sum(mu_sorted(1:p_20_ind))
        dmarg_tax_ind_q(2) = sum(dmarg_tax_sorted(p_20_ind+1:p_40_ind)*mu_sorted(p_20_ind+1:p_40_ind))/sum(mu_sorted(p_20_ind+1:p_40_ind))
        dmarg_tax_ind_q(3) = sum(dmarg_tax_sorted(p_40_ind+1:p_60_ind)*mu_sorted(p_40_ind+1:p_60_ind))/sum(mu_sorted(p_40_ind+1:p_60_ind))
        dmarg_tax_ind_q(4) = sum(dmarg_tax_sorted(p_60_ind+1:p_80_ind)*mu_sorted(p_60_ind+1:p_80_ind))/sum(mu_sorted(p_60_ind+1:p_80_ind))
        dmarg_tax_ind_q(5) = sum(dmarg_tax_sorted(p_80_ind+1:Nstate)*mu_sorted(p_80_ind+1:Nstate))/sum(mu_sorted(p_80_ind+1:Nstate))

        ! compute average transfer rates by quintile (CBO concept: total transfers of quintile divided by total income of quintile)
        avg_transfers_q(1) = sum(transfers_sorted(1:p_20_ind)*mu_sorted(1:p_20_ind))/sum(ytot_sorted(1:p_20_ind)*mu_sorted(1:p_20_ind))
        avg_transfers_q(2) = sum(transfers_sorted(p_20_ind+1:p_40_ind)*mu_sorted(p_20_ind+1:p_40_ind))/sum(ytot_sorted(p_20_ind+1:p_40_ind)*mu_sorted(p_20_ind+1:p_40_ind))
        avg_transfers_q(3) = sum(transfers_sorted(p_40_ind+1:p_60_ind)*mu_sorted(p_40_ind+1:p_60_ind))/sum(ytot_sorted(p_40_ind+1:p_60_ind)*mu_sorted(p_40_ind+1:p_60_ind))
        avg_transfers_q(4) = sum(transfers_sorted(p_60_ind+1:p_80_ind)*mu_sorted(p_60_ind+1:p_80_ind))/sum(ytot_sorted(p_60_ind+1:p_80_ind)*mu_sorted(p_60_ind+1:p_80_ind))
        avg_transfers_q(5) = sum(transfers_sorted(p_80_ind+1:Nstate)*mu_sorted(p_80_ind+1:Nstate))/sum(ytot_sorted(p_80_ind+1:Nstate)*mu_sorted(p_80_ind+1:Nstate))

        print *, ''
        print *, '**'
        print *, '** Average rates,  CBO, total system = ', avg_tax_q*100D0
        print *, '** Average rates,  ind, total system = ', avg_tax_ind_q*100D0
        print *, '** Marginal rates, ind, total system = ', marg_tax_ind_q*100D0
        print *, '** DMarginal rates, ind, total system = ', dmarg_tax_ind_q*100D0
        print *, '**'
    
    end subroutine income_distribution


    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! Distribution of tax-and-transfer rates for plots
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    subroutine implied_rates(index_ss,ngrid)

        implicit none
        
        ! Local variable declarations
        integer, intent(in) :: index_ss                 ! =1 if initial, 2 if LR
        integer, intent(in) :: ngrid                    ! number of points on which rates are computed
        real(8) :: yl_min, yl_max                       ! bounds for labor income grid
        real(8), dimension(ngrid) :: yl, yl_scaled      ! labor income grid, labor income grid relative to ymean
        real(8), dimension(ngrid) :: avg_tax_transfer   ! average tax-and-transfer rate
        real(8), dimension(ngrid) :: avg_tax            ! average tax rate
        real(8), dimension(ngrid) :: avg_transfer       ! average transfer
        real(8), dimension(ngrid) :: marg_tax_transfer  ! marginal tax-and-transfer rate
        real(8), dimension(ngrid) :: marg_tax           ! marginal tax rate
        real(8), dimension(ngrid) :: marg_transfer      ! marginal transfer rate
        integer :: iter                                 ! counter
        real(8) :: yk_mean                              ! mean capital income
        
        ! Generate labor income grid
        yl_min = 0.001d0                        ! minimum for labor income grid
        yl_max = 5d0                            ! maximum for labor income grid
        yl = LINSPACE_FET(yl_min,yl_max,ngrid)  ! labor income grid
        yl_scaled = yl*ymean                    ! rescale labor income grid
        
        open(101,  file = 'Output/yl_grid.txt',      status = 'unknown')
        write(101, *) yl
        close(101)
        
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        ! Case 1: zero capital income
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        
        ! Total tax-and-transfer rates
        do iter = 1, ngrid
            avg_tax_transfer(iter) = TAX_sc(yl_scaled(iter),0d0)/yl_scaled(iter)
            marg_tax_transfer(iter) = dTAX_sc(yl_scaled(iter),0d0)
        enddo
        
        
        if (index_ss == 1) then
            
            open(102,  file = 'Output/avg_tax_transfer_zerok_init.txt',      status = 'unknown')
            write(102, *) avg_tax_transfer
            close(102)
            
            open(103,  file = 'Output/marg_tax_transfer_zerok_init.txt',      status = 'unknown')
            write(103, *) marg_tax_transfer
            close(103)
            
            
        else
            
            open(102,  file = 'Output/avg_tax_transfer_zerok_lr.txt',      status = 'unknown')
            write(102, *) avg_tax_transfer
            close(102)
            
            open(103,  file = 'Output/marg_tax_transfer_zerok_lr.txt',      status = 'unknown')
            write(103, *) marg_tax_transfer
            close(103)
            
        endif
        
        
    end subroutine implied_rates

 

    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! Distribution hours
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    subroutine hours_distribution

        implicit none

        ! Local variable declarations
        real(8), dimension(NS) :: yl, ytot, yh                                                  ! labor income, total income, hourly wage
        real(8), dimension(NS) :: ytot_sorted, yl_sorted, yh_sorted     ! total income, labor income, tax paid, transfers (all sorted)
        real(8), dimension(NS) :: mu_sorted, mu_cumu, a_sorted     ! measure sorted, cum. measure, assets sorted, cap.+lab. tax, also sorted
        real(8), dimension(NS) :: h_sorted, h_normalized
        integer :: ind_sort(NS)                                                                 ! index for sorting
        integer :: Nstate                                                                       ! auxiliary variable equal to NS
        real(8), dimension(NS) :: muu, S1u, S2u                                                 ! auxiliary variables equal to measure, asset states, productivity states
        real(8) :: avg_h_dec(10), avg_h_q(5)
        integer :: p_20_ind, p_40_ind, p_60_ind, p_80_ind, p_90_ind                              ! indices for percentiles
        integer :: p_10_ind, p_30_ind, p_50_ind, p_70_ind                                      ! indices for percentiles
        integer :: iter, is                                                                     ! counter

        real(8) :: log_h(NS)       ! log of today's income for every state (#NS) and possible productivity tomorrow (#Nx)
        real(8) :: log_h_var          ! variance of log income
        real(8) :: log_h_sd           ! standard deviation of log income
        real(8) :: m1_log_h, m2_log_h   ! moments

        real(8) :: log_c(NS)            ! log of today's income for every state (#NS) and possible productivity tomorrow (#Nx)
        real(8) :: log_c_var            ! variance of log income
        real(8) :: log_c_sd             ! standard deviation of log income
        real(8) :: m1_log_c, m2_log_c   ! moments


        ! auxiliary variables
        Nstate = NS     ! number of individual states
        muu = mu        ! measure over individual states
        S1u = S(:,1)    ! assets for individual states
        S2u = S(:,2)    ! productivities for individual states

        ! Income variables
        yh = wge*S2u        ! labor income

        ! Normalize hours by median hours
        h_sorted = h_pol
        ind_sort = 0
        call hpsort_eps_epw(Nstate, h_sorted, ind_sort, 1D-16)
        ! sort measure and taxes paid arrays correspondingly
        mu_sorted = muu(ind_sort)                                   ! measure sorted by total income

        ! compute cumulative measure
        mu_cumu = 0d0
        mu_cumu(1) = mu_sorted(1)
        do iter = 2, Nstate
            mu_cumu(iter) = mu_cumu(iter-1) + mu_sorted(iter)
        end do

        p_10_ind = minloc(abs(mu_cumu-0.1d0), dim=1)
        p_20_ind = minloc(abs(mu_cumu-0.2d0), dim=1)
        p_30_ind = minloc(abs(mu_cumu-0.3d0), dim=1)
        p_40_ind = minloc(abs(mu_cumu-0.4d0), dim=1)
        p_50_ind = minloc(abs(mu_cumu-0.5d0), dim=1)
        p_60_ind = minloc(abs(mu_cumu-0.6d0), dim=1)
        p_70_ind = minloc(abs(mu_cumu-0.7d0), dim=1)
        p_80_ind = minloc(abs(mu_cumu-0.8d0), dim=1)
        p_90_ind = minloc(abs(mu_cumu-0.9d0), dim=1)
        !print *, 'Top quintile cutoff = ', ytot_sorted(p_80_ind)/ymean

        h_normalized = log(h_pol) - sum(mu*log(h_pol))!/Hagg!h_sorted(p_50_ind)




        ! sort hourly income array
        yh_sorted = yh
        ind_sort = 0
        call hpsort_eps_epw(Nstate, yh_sorted, ind_sort, 1D-16)
        ! sort measure and taxes paid arrays correspondingly
        mu_sorted = muu(ind_sort)                                   ! measure sorted by total income
        h_sorted = h_normalized(ind_sort)

        ! compute cumulative measure
        mu_cumu = 0d0
        mu_cumu(1) = mu_sorted(1)
        do iter = 2, Nstate
            mu_cumu(iter) = mu_cumu(iter-1) + mu_sorted(iter)
        end do

        p_10_ind = minloc(abs(mu_cumu-0.1d0), dim=1)
        p_20_ind = minloc(abs(mu_cumu-0.2d0), dim=1)
        p_30_ind = minloc(abs(mu_cumu-0.3d0), dim=1)
        p_40_ind = minloc(abs(mu_cumu-0.4d0), dim=1)
        p_50_ind = minloc(abs(mu_cumu-0.5d0), dim=1)
        p_60_ind = minloc(abs(mu_cumu-0.6d0), dim=1)
        p_70_ind = minloc(abs(mu_cumu-0.7d0), dim=1)
        p_80_ind = minloc(abs(mu_cumu-0.8d0), dim=1)
        p_90_ind = minloc(abs(mu_cumu-0.9d0), dim=1)
        !print *, 'Top quintile cutoff = ', ytot_sorted(p_80_ind)/ymean

        ! compute average hours  by decile
        avg_h_dec(1) = sum((h_sorted(1:p_10_ind))*mu_sorted(1:p_10_ind))/sum(mu_sorted(1:p_10_ind))
        avg_h_dec(2) = sum((h_sorted(p_10_ind+1:p_20_ind))*mu_sorted(p_10_ind+1:p_20_ind))/sum(mu_sorted(p_10_ind+1:p_20_ind))
        avg_h_dec(3) = sum((h_sorted(p_20_ind+1:p_30_ind))*mu_sorted(p_20_ind+1:p_30_ind))/sum(mu_sorted(p_20_ind+1:p_30_ind))
        avg_h_dec(4) = sum((h_sorted(p_30_ind+1:p_40_ind))*mu_sorted(p_30_ind+1:p_40_ind))/sum(mu_sorted(p_30_ind+1:p_40_ind))
        avg_h_dec(5) = sum((h_sorted(p_40_ind+1:p_50_ind))*mu_sorted(p_40_ind+1:p_50_ind))/sum(mu_sorted(p_40_ind+1:p_50_ind))
        avg_h_dec(6) = sum((h_sorted(p_50_ind+1:p_60_ind))*mu_sorted(p_50_ind+1:p_60_ind))/sum(mu_sorted(p_50_ind+1:p_60_ind))
        avg_h_dec(7) = sum((h_sorted(p_60_ind+1:p_70_ind))*mu_sorted(p_60_ind+1:p_70_ind))/sum(mu_sorted(p_60_ind+1:p_70_ind))
        avg_h_dec(8) = sum((h_sorted(p_70_ind+1:p_80_ind))*mu_sorted(p_70_ind+1:p_80_ind))/sum(mu_sorted(p_70_ind+1:p_80_ind))
        avg_h_dec(9) = sum((h_sorted(p_80_ind+1:p_90_ind))*mu_sorted(p_80_ind+1:p_90_ind))/sum(mu_sorted(p_80_ind+1:p_90_ind))
        avg_h_dec(10) = sum((h_sorted(p_90_ind+1:Nstate))*mu_sorted(p_90_ind+1:Nstate))/sum(mu_sorted(p_90_ind+1:Nstate))

        ! By quintile
        ! compute average hours  by decile
        avg_h_q(1) = sum((h_sorted(1:p_20_ind))*mu_sorted(1:p_20_ind))/sum(mu_sorted(1:p_20_ind))
        avg_h_q(2) = sum((h_sorted(p_20_ind+1:p_40_ind))*mu_sorted(p_20_ind+1:p_40_ind))/sum(mu_sorted(p_20_ind+1:p_40_ind))
        avg_h_q(3) = sum((h_sorted(p_40_ind+1:p_60_ind))*mu_sorted(p_40_ind+1:p_60_ind))/sum(mu_sorted(p_40_ind+1:p_60_ind))
        avg_h_q(4) = sum((h_sorted(p_60_ind+1:p_80_ind))*mu_sorted(p_60_ind+1:p_80_ind))/sum(mu_sorted(p_60_ind+1:p_80_ind))
        avg_h_q(5) = sum((h_sorted(p_80_ind+1:Nstate))*mu_sorted(p_80_ind+1:Nstate))/sum(mu_sorted(p_80_ind+1:Nstate))

        print *, ''
        print *, 'log(Hours/mean hours) by wage'
        !print *, '** by decile ', avg_h_dec
        print *, '** by quintile ', avg_h_q


        ! sort by hours
        h_sorted = h_normalized
        ind_sort = 0
        call hpsort_eps_epw(Nstate, h_sorted, ind_sort, 1D-16)
        ! sort measure and taxes paid arrays correspondingly
        mu_sorted = muu(ind_sort)                                   ! measure sorted by total income

        ! compute cumulative measure
        mu_cumu = 0d0
        mu_cumu(1) = mu_sorted(1)
        do iter = 2, Nstate
            mu_cumu(iter) = mu_cumu(iter-1) + mu_sorted(iter)
        end do

        p_10_ind = minloc(abs(mu_cumu-0.1d0), dim=1)
        p_20_ind = minloc(abs(mu_cumu-0.2d0), dim=1)
        p_30_ind = minloc(abs(mu_cumu-0.3d0), dim=1)
        p_40_ind = minloc(abs(mu_cumu-0.4d0), dim=1)
        p_50_ind = minloc(abs(mu_cumu-0.5d0), dim=1)
        p_60_ind = minloc(abs(mu_cumu-0.6d0), dim=1)
        p_70_ind = minloc(abs(mu_cumu-0.7d0), dim=1)
        p_80_ind = minloc(abs(mu_cumu-0.8d0), dim=1)
        p_90_ind = minloc(abs(mu_cumu-0.9d0), dim=1)
        !print *, 'Top quintile cutoff = ', ytot_sorted(p_80_ind)/ymean

        ! compute average hours  by decile
        avg_h_dec(1) = sum((h_sorted(1:p_10_ind))*mu_sorted(1:p_10_ind))/sum(mu_sorted(1:p_10_ind))
        avg_h_dec(2) = sum((h_sorted(p_10_ind+1:p_20_ind))*mu_sorted(p_10_ind+1:p_20_ind))/sum(mu_sorted(p_10_ind+1:p_20_ind))
        avg_h_dec(3) = sum((h_sorted(p_20_ind+1:p_30_ind))*mu_sorted(p_20_ind+1:p_30_ind))/sum(mu_sorted(p_20_ind+1:p_30_ind))
        avg_h_dec(4) = sum((h_sorted(p_30_ind+1:p_40_ind))*mu_sorted(p_30_ind+1:p_40_ind))/sum(mu_sorted(p_30_ind+1:p_40_ind))
        avg_h_dec(5) = sum((h_sorted(p_40_ind+1:p_50_ind))*mu_sorted(p_40_ind+1:p_50_ind))/sum(mu_sorted(p_40_ind+1:p_50_ind))
        avg_h_dec(6) = sum((h_sorted(p_50_ind+1:p_60_ind))*mu_sorted(p_50_ind+1:p_60_ind))/sum(mu_sorted(p_50_ind+1:p_60_ind))
        avg_h_dec(7) = sum((h_sorted(p_60_ind+1:p_70_ind))*mu_sorted(p_60_ind+1:p_70_ind))/sum(mu_sorted(p_60_ind+1:p_70_ind))
        avg_h_dec(8) = sum((h_sorted(p_70_ind+1:p_80_ind))*mu_sorted(p_70_ind+1:p_80_ind))/sum(mu_sorted(p_70_ind+1:p_80_ind))
        avg_h_dec(9) = sum((h_sorted(p_80_ind+1:p_90_ind))*mu_sorted(p_80_ind+1:p_90_ind))/sum(mu_sorted(p_80_ind+1:p_90_ind))
        avg_h_dec(10) = sum((h_sorted(p_90_ind+1:Nstate))*mu_sorted(p_90_ind+1:Nstate))/sum(mu_sorted(p_90_ind+1:Nstate))

        ! By quintile
        ! compute average hours  by decile
        avg_h_q(1) = sum((h_sorted(1:p_20_ind))*mu_sorted(1:p_20_ind))/sum(mu_sorted(1:p_20_ind))
        avg_h_q(2) = sum((h_sorted(p_20_ind+1:p_40_ind))*mu_sorted(p_20_ind+1:p_40_ind))/sum(mu_sorted(p_20_ind+1:p_40_ind))
        avg_h_q(3) = sum((h_sorted(p_40_ind+1:p_60_ind))*mu_sorted(p_40_ind+1:p_60_ind))/sum(mu_sorted(p_40_ind+1:p_60_ind))
        avg_h_q(4) = sum((h_sorted(p_60_ind+1:p_80_ind))*mu_sorted(p_60_ind+1:p_80_ind))/sum(mu_sorted(p_60_ind+1:p_80_ind))
        avg_h_q(5) = sum((h_sorted(p_80_ind+1:Nstate))*mu_sorted(p_80_ind+1:Nstate))/sum(mu_sorted(p_80_ind+1:Nstate))

        print *, ''
        print *, 'log(Hours/mean hours) by hours'
        !print *, '** by decile ', avg_h_dec
        print *, '** by quintile ', avg_h_q



        log_h = log(h_pol)
        m1_log_h = sum(log_h*muu)
        m2_log_h = sum(muu*(log_h-m1_log_h)**2d0)
        log_h_var = m2_log_h
        log_h_sd = sqrt(m2_log_h)
        print *, ''
        print *, "Variance of log h = ", log_h_var


        log_c = log(c_pol)
        m1_log_c = sum(log_c*muu)
        m2_log_c = sum(muu*(log_c-m1_log_c)**2d0)
        log_c_var = m2_log_c
        log_c_sd = sqrt(m2_log_c)
        print *, ''
        print *, "Variance of log c = ", log_c_var


    end subroutine hours_distribution

   
end module
